Ab Initio Potentials for the Ground S0 and the First Electronically Excited Singlet S1 States of Benzene–Helium with Application to Tunneling Intermolecular Vibrational States

We present new ab initio intermolecular potential energy surfaces for the benzene–helium complex in its ground (S0) and first excited (S1) states. The coupled-cluster level of theory with single, double, and perturbative triple excitations, CCSD(T), was used to calculate the ground state potential. The excited state potential was obtained by adding the excitation energies S0 → S1 of the complex, calculated using the equation of motion approach EOM-CCSD, to the ground state potential interaction energies. Analytical potentials are constructed and applied to study the structural and vibrational dynamics of benzene–helium. The binding energies and equilibrium distances of the ground and excited states were found to be 89 cm–1, 3.14 Å and 77 cm–1, 3.20 Å, respectively. The calculated vibrational energy levels exhibit tunneling of He through the benzene plane and are in reasonable agreement with recently reported experimental values for both the ground and excited states [Hayashi, M.; Ohshima, Y. J. Phys. Chem. Lett.2020, 11, 9745]. Prospects for the theoretical study of complexes with large aromatic molecules and He are also discussed.


I. INTRODUCTION
The molecular complexes of aromatic molecules with He atoms (ArmHe) stand out from other complexes of aromatic molecules with heavier Rg atoms (ArmRg) because they are extremely weakly bound.The potential energy depth for ArmHe is about three times less than for an analogous ArmNe complex and about eight times less than for ArmAr.−3 The zero-point energies of the complexes with He are lying above half of the binding energies of He; in this case, the wave functions (WF) of He motion are spread through the monomer plane, and the vibrational levels exhibit tunneling splitting.−15 Recently, Hayashi and Ohshima reported very high-resolution UV excitation spectra of BzHe in which the splitting of the two lowest vibrational excitations for the excited state S 1 was detected. 16Some small splitting in the rotational transitions was also noticed earlier for the pyridine−He complex, 17 which allegedly stems from He delocalization through the monomer plane.
The ArmHe complexes have significant implications in the context of spectroscopy of aromatic molecules embedded in helium droplets.This has been extensively studied in the literature, including refs 18−22.Electronic spectroscopy of molecules in superfluid helium nanodroplets reveals interesting effects, such as vibronic bands of 0 0 0 zero phonon lines (ZPL) followed by broader "phonon wings" (PW) on the blue side, as described in refs 23 and 24.The existence of superfluidity in helium nanodroplets was first reported in ref 23 based on the gap between these bands.The electronic transition S 0 → S 1 of the aromatic molecules to the first excited single state has been observed in helium droplets for a series of aromatic molecules, as reported in refs 14, 19, and 25−31.Furthermore, a study of the impact of the rotation of a benzene molecule on the helium density distribution in a helium droplet was conducted in ref 32.−41 In some cases, these potentials were adjusted to reproduce experimental data, but the accuracy of the potentials should be confirmed by more reliable PESs based on ab initio calculations.It is possible that such potentials may differ significantly from ab initio potentials and inaccurately reproduce experimental results like it takes place for BzKr. 42he ab initio ground state PES for BzHe was previously reported, 5 aiming to complement the experimental data on rotational constants and Raman scattering coefficients also presented in that study.The reported potential was then used to study the vibrational structure of the complex, and its accuracy was recently confirmed by experiment. 16However, our tests of the reported potential revealed some unphysical behavior in the long-range region of the potential.−45 Therefore, a new potential for the BzHe complex in its ground state with correct behavior in the long-range region is needed.Previous theoretical studies for other ArmHe complexes where an initio PES was reported are rather scarce, and no vibrational structure was studied. 46,47lthough most experimental studies of ArmRg complexes focus on the excited state S 1 , the ground state potential can still be useful for assigning the experimental spectrum due to the similar shapes of the potentials for both states.Nevertheless, theoretical ab initio potentials for the excited state S 1 have been reported for a series of ArmAr complexes in previous studies, 48−51 which unsurprisingly gave better agreement with experimental vibrational energies observed for the S 1 state than their ground state counterparts.In these studies, the excited state E int was calculated by adding the adiabatic excitation energy found with the single-reference linear-response coupled cluster (LR-CCSD) method 52,53 to the ground state CCSD(T) potential.It is worth noting that the equation-of-motion coupled cluster (EOM-CCSD) method, 54 which is another way of treating electronically excited states, is equivalent to LR-CCSD when it comes to calculating excitation energy values.
The objective of this work is to construct a new potential for the BzHe ground state with the correct asymptotic behavior and at a higher level of theory.Additionally, the potential for the excited state S 1 will be reported.Then, both potentials will be employed to study the tunneling motion through the monomer plane using the appropriate theoretical approach.We expect that our new potentials will be useful for the community.
The structure of our article is as follows: First, the electronic structure methods used in this work are described in Section II.Next, the analytical form of the potential is presented in Section III, followed by an analysis of the vibrational structure in Section IV.Finally, future prospects and challenges in this field are discussed in Section V.

II. AB INITIO METHODS
The intermolecular Cartesian coordinates (x, y, z) describing the motion of He are chosen as follows.The reference point is at the benzene center of mass, the z-axis coincides with the C 6 symmetry axis of benzene, and the y-axis crosses the midpoint of two C−C bonds.The topology of potential PES is typical for BzRg complexes; it has two symmetrical global minima M z above and below the benzene plane and six symmetrical local minima M y in the monomer plane.The minimum-energy paths (MEPs) connect the minima through two kinds of saddle points, S x and S yz , six of each type.These coordinates can also be transformed into the standard spherical coordinates (r, θ, φ).
For the calculation of the ground state S 0 potential energy E int , we used the conventional CCSD(T) method with the frozen core approximation (FC) and the counterpoise correction, 55 along with Dunning's basis sets aug-cc-pVXZ (X = D, T, Q). 56,57 The basis sets were augmented with 3s3p2d1f1g midbond functions, 58,59 with exponents of 0.9, 0.3, and 0.1 for s and p functions, 0.6 and 0.2 for f functions, 0.3 for the f function, and 0.25 for the g function.The notation aug-cc-pVXZ+ was used for these basis sets with midbonds.The convergence of E int was investigated near the critical points M z , M y , and S x by local interpolation, and it was found that the difference between augcc-pVTZ+ and aug-cc-pVQZ+ in the coordinates and energy values was very small, as shown in Table 1.Therefore, the CCSD(T)/aug-cc-pVTZ+ level of theory was used for the calculation of E int over the entire physically meaningful range of the intermolecular coordinates.It should be noted that the authors of ref 5 employed CCSD(T)/aug-cc-pVDZ+ and did not calculate some important single-point energy values in the long-range region, which led to unphysical behavior of the potential.Details of the unphysical behavior of their potential are given in the Supporting Information.
The excited state S 1 of the BzHe complex dissociates to monomers with the following reaction: The excited state intermolecular interaction energy E int is calculated using the supermolecular approach: The excited state energies E BzHe S 1 and E Bz S 1 can be calculated by the EOM-CCSD method.However, E int found with CCSD for the ground state and with EOM-CCSD for the excited state is sensitive to the absence of important triple contributions.As shown in Table 1, the binding energy D e for the ground state potential is severely underestimated when compared to the The Journal of Physical Chemistry A CCSD(T) values.On the other hand, employing the more accurate EOM-CCSD(T) 60 or EOM-CCSDT 61 in eq 2 would be too time-consuming in a reasonable basis set for this complex.−51 In this approach, the adiabatic excitation energies for BzHe and Bz are found in the dimer basis as the difference between EOM-CCSD (or LR-CCSD) energies for S 1 calculated using the excited state geometry of benzene and the ground state CCSD energies for S 0 calculated using the ground state benzene geometry.Then, the difference between the excitation energies for BzHe and Bz is added to the ground state CCSD(T) potential, giving the excited state E int .This value has to be corrected by the difference between He energies found in the dimer basis set in the excited and ground state benzene geometries to eliminate the basis set superposition error (BSSE).It is convenient to write this approach as follows: The method described by eq 3 will be denoted as EOM-CCSD(T)* in this work.To calculate the excited state correction, we tested two basis sets, namely aug-cc-pVDZ and aug-cc-pVTZ.As shown in Table 1, the difference in energy values calculated with these two basis sets is almost negligible.Therefore, we used the smaller basis set and a C 1 point group symmetry at the same He positions as for the ground state calculations.The overall number of He positions was 349 in all the calculations.The number is excessive because, after preliminary calculations, an additional set of He positions was added in the region above the zero-point vibrational energy to ensure an accurate description of the vibrational energy levels.
The binding energy D e can be found as −V(0,0,z e ), where z e is the equilibrium distance.This value in the excited S 1 state is 12 cm −1 lower than its ground state counterpart, whereas the BzAr complex in the excited S 1 state is characterized by a larger value of D e than in the ground state. 48,49A comparison of the selected 1-dimensional cuts of the ground and excited state potentials is shown in Figure 1.As one can notice, the excited state potential is flatter in the vicinity of the global minimum and characterized by a lower value of the potential barrier at the saddle point S yz .In

The Journal of Physical Chemistry A
the monomer plane, the E int for the ground state is slightly flatter than the excited state one.
We used the same ground state benzene geometry from ref 62 −65 The ground state bond lengths were 1.3915 Å for CC and 1.080 Å for CH bonds.The excited state S 1 geometry was taken from ref 66 with bond lengths of 1.430 Å for CC and 1.085 Å for CH bonds.All electronic structure calculations were performed using the Molpro software package. 67The Supporting Information contains the calculated energy values.

III. ANALYTICAL PES
There are a few possible ways to convert the calculated set of single-point E int to an analytical form.One way to use the expansion of E int in spherical harmonics basis, which in the case of a nonlinear molecule and atoms has the following form: where m = −l, −l + 1, ..., l and R = (R,θ,φ) is the vector connecting the center of mass of a nonlinear molecule and an atom.In eq 4, Ω lm (θ,φ) stands for the normalized tesseral (real) harmonics basis. 68In the case when the nonlinear molecule has It is convenient to divide the functions χ lm (R) into short-and long-range parts χ lm (R) = χ lm sr (R) + χ lm lr (R).The functions χ lm sr (R) and χ lm lr (R) must vanish in the long and short ranges, respectively.In this case the potential can be written as follows: The long-range potential can be presented using the form where it is reasonable to take 1/R n up to n = 10 for accurate reproduction of the long-range potential and the terms with (l, m) = (0, 0), (2, 0), (4, 0), (6, 0), (6, 6) due to benzene symmetry.The term with n = 6 describes the dispersion coupling of induced dipole polarizabilities of the monomers and the terms with n = 8 dipole−quadrupole and n = 10 quadrupole− quadrupole couplings. 65In order to guarantee vanishing of V lr (R) in the short-range regions, the expansion terms have to multiplied by damping functions f n (R).A popular choice is the Tang−Toennies functions 69 defined as . We employ another type of damping function as explained below.
When it comes to the short-range interaction potential, we found that it is more convenient to use a many-body form of potential with variables r k and functions w(r k ) = 1 − e −a k (r k −r ke )/r ke .Here r k is the distance between an atom and kth atom of an aromatic molecule, and a k and r ke are some parameters found from nonlinear optimization fitting algorithm.−65 It can be written as follows: In eq 7 the two-body terms V 2 are defined as and three-body terms as The constants V 0 and W 0 are chosen manually to reproduce ab initio E int with maximum accuracy.
For such form of the analytical potential as given in eq 7, we used another kind of damping functions for the short-and long-range parts, denoted as h sr (R) and h lr (R), used in eq 6 instead of f n lr (R): First, we found the set of the long-range coefficients C n l,m by fitting E int to eq 6 for R ⩾ 8 Å for both the ground and excited states PESs.V lr (R) reproduces the long-range E int accurately with root-mean-square error (RMSE) of order of 0.001 cm −1 .The obtained C n l,m agree well with their counterparts found from the first-principles using monomer properties such as multipole moments and static/dynamics polarizabilities. 70More details on the comparison of C n l,m obtained by different methods can be found in the Supporting Information.
Then, the long range part was kept frozen, and the parameters {c i , c ij , a k , r ke , R 0 , γ} of V sr (R) in eq 7 were found using the nonlinear Levenberg−Marquardt algorithm. 71,72The final potential described by eq 5 is characterized by RMSE equal to 0.2 (0.6) cm −1 in the energy range (−D e , 0) and 1.0 (0.6) cm −1 in the range (0, 500 cm −1 ) for S 0 (S 1 ) PES.Both potentials have correct asymptotic behavior in the short-and long-range regions.
Unlike BzRg complexes with heavier Rg atoms, the BzHe potential is more anisotropic with respect to φ variable in the vicinity of the M z point.The saddle points S yz are located at a larger distance from the monomer plane than M z .The potential barrier at S yz is only 58.9 cm −1 for the S 1 PES that allows tunneling motion of He.As one can see in Figure 2a−d, the delocalization channels connect the minima M z above and below the monomer plane.The saddle points S x and S yz are connected with the minima M z and M y through the minimum-energy path (MEP), shown only for the S x minima in Figure 2.
A Fortran 90 program for the calculation of E int for both PESs at a given position of He is given in the Supporting Information.

IV. TUNNELING VIBRATIONAL STATES
When a heavier rare gas atom is used in a complex with benzene, the tunneling motion of Rg can often be neglected for the lowest vibrational states, allowing for an approach using Cartesian coordinates to solve the nuclear 3D vibrational Schrodinger equation, 42,63 in which the vibrational basis set functions are located at one global minimum M z .In principle, it is possible to use the approach to study the tunneling effect by locating the basis set functions either on both minima M z or at the center of the coordinate frame.Instead, we employed an alternative approach that utilizes spherical coordinates (r, θ, φ) proposed in refs 73 and 74.Details of the numerical procedure are described in ref 63.Within the procedure, the vibrational WF are represented as a series of product basis functions constructed The Journal of Physical Chemistry A from stretching−bending and free internal rotation basis functions.Table 2 contains the lowest vibrational energy levels for BzHe in both considered electronic states.The results in Table 2 are given with labels for the symmetry group G 24 used in ref 5 for compatibility of the results.The zero-point energies are given relative to dissociation and are found as the difference between the D e values and the ground intermolecular vibrational state.Figure 3 contains cuts of the vibrational wave functions for the selected lowest states at φ = 0.The redistributions below and above the monomer plane clearly indicate the presence of quantum tunneling effect.
In contrast to BzRg with heavier than He atoms, 64,65 the assignment of the vibrational states in terms of the harmonic oscillator quantum numbers (n s , n b , l) for the singly degenerate stretching mode n s and the doubly degenerate bending mode (n b , l) can only be done formally.The states 2A 1 ′/2A 1 ″ and 1E 1 ′/ 1E 1 ″ in Table 2 are approximately identified as the lowest stretching and bending modes, respectively.
The results for the S 0 PES are very close to those obtained by Lee et al., though not entirely identical.Several reasons account for this disparity.First, we utilized a larger basis set for our CCSD(T) calculations, encompassing a wider set of He positions.Second, our potential exhibits correct asymptotic behavior, which may affect the higher vibrational states close to the dissociation limit.Third, the theoretical approach to treat the vibrational problem differs from that employed by Lee et al.
Our results for the S 1 PES revealed the following: the value of the lowest vibrational state, corresponding to the nominal bending mode, is almost 2 cm −1 closer to the experimental value than its S 0 counterpart.The splitting value of the vibrational energy levels extracted from the S 1 PES is in good agreement with the experimental estimation of <0.02.However, the subsequent state, corresponding to nominal stretching for the S 1 state, exhibit less favorable agreement with the experimental values of the vibrational energy.
We have also calculated the values of averaged intermolecular distances R 0 for the zero-point vibrational energies, which are 3.590 and 3.604 Å in the S 0 and S 1 states, respectively.These values agree reasonably well with their counterparts derived from the experiment: 3.602 and 3.665 Å, respectively. 6he vibrational energy levels for the C 6 D 6 isotopomer for both electronic states, as well as the energy levels for higher vibrational energies than those given in Table 2, are reported in the Supporting Information.

V. CONCLUSIONS
We have presented a study of the BzHe complex in its ground S 0 and the first singlet excited electronic states S 1 .The conventional The Journal of Physical Chemistry A CCSD(T) method with augmented Dunning's triple-ζ basis set supplemented by midbond functions has been employed to calculate the ground state E int .We obtain the excited state E int by adding the adiabatic correction, found using the EOM-CCSD method, to the ground state CCSD(T) potential.A nonlinear optimization algorithm is used to construct analytical potentials for each state.Our potential for the BzHe ground state has the correct behavior in the long-range region compared to the previously reported potential. 5oth potentials exhibit similar shapes, although the excited state potential appears visibly flatter in the vicinity of the global minimum.We apply these potentials to solve the 3D vibrational Schrodinger equation, which describes the motion of an atom with respect to a plane symmetric top molecule using spherical coordinates.This approach enables us to study the tunneling motion of an atom through the monomer plane.The calculated vibrational energy levels exhibit splitting, and the vibrational wave functions are dispersed throughout the monomer plane.The vibrational energy values obtained from these potentials reasonably agree with recently observed experimental values 16 for both electronic states, S 0 and S 1 .However, the somewhat surprising disagreement in the vibrational energy value corresponding to nominal stretching for the S 1 state necessitates further investigation.This discrepancy may stem from inaccuracies in the single-reference excited states methods, such as EOM-CCSD, particularly in regions closer to the dissociation limit, and multireference methods may help to remove this inaccuracy. 75 problem of great importance is the calculation of PES for complexes containing large aromatic molecules such as tetracene and coronene with He. 22 Unfortunately, the "goldstandard" CCSD(T) becomes too time-consuming for such complexes.Therefore, there is a need for cheaper methods that provide similar accuracy at much less computational cost.We recently showed that symmetry-adapted perturbation theory based on density functional description of monomer properties (DFT-SAPT) can be reduced to an approximate model 76 in which the time-consuming dispersion term can be replaced by its empirical counterpart. 77Another alternative is the secondorder Møller−Plesset perturbation theory which, according to our preliminary tests, gives good agreement of E int with CCSD(T) for ArmHe and ArmNe complexes thanks to fortuitous error cancellation.A detailed study of the potentials for complexes with larger aromatic molecules will be presented in a separate study.
Figure S1 and Tables S1

Notes
The author declares no competing financial interest.

■ ACKNOWLEDGMENTS
The author thanks the PL-GRID project for granting access to their computational resources and the PL-GRID team for helpful software support.Professor Jan Makarewicz is acknowledged for stimulating the author's interest in the studied complex.

Figure 1 .
Figure 1.1D cuts of the ground and excited state potentials along the x-, y-, and z-axes and the line x = 0, z = z e .The solid lines denote the potential cuts normalized to zero at infinity.ΔE BzHe and ΔE Bz are zero−zero adiabatic excitation energies for BzHe and Bz, respectively.

Figure 2 .
Figure 2. PES cuts for BzHe complex in the xy plane at z = z e (a), at z = 0 (b), at z = 0.5z e (c), and at z = 1.1z e (d).The coordinates are in Å, and the contour values are in cm −1 .

Figure 3 .
Figure 3. 2D cuts of the normalized WFs of BzHe for the lowest vibrational states for the S 0 potential at φ = 0.The WF isovalues differ by 0.2, where the maximum value of 1.0 corresponds to the WF maximum.

Table 1 .
Values (α, E int ), Where α Is the z, y, and x Coordinates (in Å) of the Critical Points M z , M y , and S x , Respectively, and (y, z, E int ) for the Saddle Point S yz , and E int Is the Energy Value at These Points (in cm −1 ), Calculated Using Different Methods for the Ground S 0 and the Excited S 1 States of BzHe.Notations for the Basis Sets avdz=aug-cc-pVDZ and avtz=aug-cc-pVTZ Are Used.Found from analytical PES; see Section III.b Basis set used for S 0 CCSD(T)/basis set used for EOM-CCSD S 1 ; see eq 3. a